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ABSTRACT 

A proper test of Modified Newtonian Dynamics (MOND) in systems of non-trivial ge- 
ometries depends on modelling subtle differences in several versions of its postulated 
theories. This is especially important for lensing and dynamics of barely virialised 
galaxy clusters with typical gravity of scale ~ ao ~ lAs~ 2 . The original MOND for- 
] mula, the classical single field modification of the Poisson equation, and the multi-field 

\Q ■ general relativistic theory of Bekenstein (TeVeS) all lead to different predictions as we 

c"^i ' stray from spherical symmetry In this paper, we study a class of analytical MON- 

, Dian models for a system with a semi-Hernquist baryonic profile. After presenting the 

' analytical distribution function of the baryons in spherical limits, we develop orbits 

and gravitational lensing of the models in non-spherical geometries. In particular, we 
can generate a multi-centred baryonic system with a weak lensing signal resembling 
I 1 that of the merging galaxy cluster IE 0657-56 with a bullet-like light distribution. 

We finally present analytical scale- free highly non-spherical models to show the subtle 
+-> ■ differences between the single field classical MOND theory and the multi-field TeVeS 

^ ' theory. 
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1 INTRODUCTION 

While luminous models for the central parts of galaxies do not usually require a dark matter (DM) component, massive halos 
of DM must be taken into account if one wants to construct realistic potential-density pairs for an entire galaxy, in order to 
reproduce the observed nearly-flat galactic rotation curves. 

Curiously, there is a considerable body of evidence that the galactic mass profiles of baryonic and dark matter are not 
uncorrelated (McGaugh 2005). The correlation between the Newtonian gravity of the baryons and the overall gravity g 
(baryons plus DM) can be loosely parameterized by Milgrom's (1983) empirical relation 

Ks/ao)g = gN, (1) 

where the interpolating function n(x) is a function which runs smoothly from fi(x) — x at x <C 1 to jJ,(x) = 1 at i > 1 
with a dividing gravity scale ao ~ lAs -2 at the transition. It is a very serious challenge for cold dark matter simulations to 
reproduce this empirical relation. On the other hand, this relation can be interpreted as a modification of the Newton-Einstein 
gravitational law in the ultra- weak field regime, with no actual need for dark matter (at the galaxy scale). This provocative 
idea was taken as the basis for the MOND theory: however, using Eq.(l) alone to transform gN into g does not provide a 
respectable gravitational force preserving energy and angular momentum. 

Bekenstein & Milgrom (1984) then suggested modifying the Poisson Equation in order to produce a non-relativistic MOND 
potential. This proposal was recently refined by Bekenstein (2004) who presented a Lorentz-covariant theory of gravity, dubbed 
TeVeS, yielding MONDian behaviour in the appropriate limit. However, there is a subtle difference between the non-relativistic 
formulation of MOND and the relativistic TeVeS, namely the number of fields involved: in MOND, the potential is modified 
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directly in order to trigger the MOND phenomenology, while in TeVeS, a scalar field is added to the traditional Newtonian 
potential. These two descriptions are equivalent only in highly symmetric systems (spherical or cylindrical symmetry). For 
disk galaxies, Eq.£Ql is a good approximation of the MOND and TeVeS gravitational theories since the additional curl field 
(see is small when solving the modified Poisson equation for the potential or the scalar field (Brada & Milgrom 1995). 
For this reason, Eq.Q has been used with confidence to fit the rotation curves of an impressive list of external galaxies with 
remarkable accuracy (Sanders & McGaugh 2002). However, very little work has been carried out to study the actual effect 
of the curl field in the non-relativistic MOND theory (Brada & Milgrom 1995, 1999; Ciotti, Londrillo & Nipoti 2006), while 
the quantitative difference between the predictions of MOND and TeVeS has not been studied at all. The issue has become 
urgent as recent evidence against MOND is largely based on multi-centred systems such as satellites of the Milky Way (Zhao 
2005) and the bullet cluster of merging galaxies (Clowe et al. 2004). 

In this paper, after recalling some basic formulae of MOND and TeVeS (3U, we propose a set of parametric interpolating 
functions that are physical in TeVeS (3HJ. Then we present spherical analytical potential-density pairs for the baryon distri- 
bution in early type galaxies, galactic bulges and dwarf spheroidals, valid in MOND as well as in TeVeS (*SJ. For this model 
we work out the density, potential, circular velocity, isotropic and anisotropic distribution functions. We also show a simple 
analytical expression for the gravitational bending angle. Then we present the combination of such models in multi-centred 
systems for the classical MOND approximation(f|I>J; this provides an indication of the kind of result we could expect from a 
rigorous modelling of the bullet cluster of galaxies (Clowe et al. 2004). More generally, we then show how to extend those 
spherical models to oblate ones (^J. In those non-spherical situations, the density corresponding to a given potential differs in 
the classical MOND and in TeVeS. To illustrate this, we study the effects of the flattening of the potential in the special case of 
scale- free oblate one-dimensional models (JTJ, and we explicitly show the subtle differences between the different theoretical 
frameworks. Note that Milgrom (1994) has also suggested that MOND could have a modified inertia basis rather than a 
modified gravity basis; in that case the original MOND formula is correct for circular orbits in galaxies, the theory would 
become strongly non-local, the conservation laws would become unusual, and the potential-density approach used hereafter 
would not apply to that framework. 



2 MOND AND TEVES 

In the aquadratic Lagrangian theory of MOND by Bekenstein & Milgrom (1984), the Poisson equation reads 

V.[mV$] = V 2 $iv = 4tyG P , (2) 

where ^i(|V$|/ao) is the same interpolating function as in Eq.Q. We then have 

Kfl/ao) g = gN + V x H. (3) 

The value of the curl field depends on the boundary conditions, but vanishes in spherical symmetry where Gauss' theorem 
applies and Eq.Q is recovered. In realistic geometries, the curl field is non-zero but small (see Brada & Milgrom 1995, 1999; 
Ciotti et al. 2006), leading to small differences when computing the rotation curves of spiral galaxies. 

On the other hand, Bekenstein's relativistic MOND (Bekenstein 2004) is a tensor-vector-scalar (TeVeS) theory: the tensor 
is an Einstein metric g a p out of which is built the usual Einstein-Hilbert action. U a is a dynamical normalized vector field 
(g al3 U a Uf3 — —1), and <j> is a dynamical scalar field. The action is the sum of the Einstein-Hilbert action for the tensor g a g, 
the matter action, the action of the vector field U a , and the action of the scalar field <f). Einstein-like equations are obtained 
for each of these fields by varying the action w.r.t. each of them. 

In TeVeS, the physical metric near a quasi-static galaxy is given by the same metric as in General Relativity, with the 
Newtonian potential <I>jv replaced by the total potential 

$ = H$jv + 0, (4) 

where S ~ 1. This means that the scalar field <f> plays the role of the dark matter gravitational potential. The Einstein-like 
equation for the scalar field relates it to the Newtonian potential $jv (generated by the baryonic density p) through the 
equation 

V.[m s V<£] = V 2 $jv = 4nGp, (5) 
where jj, s is a function of the scalar field strength g s — |V0|, and derives from a free function in the action of the scalar field. 



3 THE INTERPOLATING FUNCTIONS 

In spherical symmetry, we have 

fi s g s = fi(g s + gN) = gN, (6) 
where p is the interpolating function of MOND, thus related to p a by 
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The standard interpolating function that has been used for twenty years to fit the rotation curves with Eq.(f) is 

H{x) = —^=. (8) 
Vl + x 

However, Zhao & Famaey (2006) have shown that this function, or rather the corresponding function fi s derived from Eq.0, 
is not physical in the framework of TeVeS; the function /j, 3 is multi- valued. For this reason, we will use hereafter a parametric 
set of interpolating functions that are physical in TeVeS: 

2r 

fi(x) = — , < a < 1, (9) 

1 + (2 - a)x + ^(1 - ax) 2 + 4a; 

where x — g/ao- The corresponding scalar field functions are 
1 — as 

where s = g s /ao- The a — case corresponds to the toy model proposed by Bekenstein (2004) in weak and intermediate 
gravity. Under the approximation of Eq.(l), the a = 1 model has been shown by Famaey & Binney (2005) and Zhao & 
Famaey (2006) to be a better fit to the rotation curves of galaxies than the a = case. The values < a < 1 have not yet 
been explored in real galaxies, and will be the subject of another paper (Famaey, Gentile & Zhao 2006, in preparation). 

These functions and fi s will lead to the same gravitational behaviour in spherical symmetry. However, note that Eq.JTJ 
is not valid in a more general geometry: the Newtonian force, the MOND force of Bekenstein & Milgrom (1984) and the TeVeS 
force are no longer parallel. The curl field obtained when solving the equation for the scalar field cf> (Eq. |SJ will be different 
than the one obtained when solving for the full $ in Eq.(|5J. This will be illustrated in §7. 



4 SPHERICAL POTENTIAL-DENSITY PAIRS IN MOND/TEVES 

In this section, we consider only the a = 1 case, known to yield excellent fits to the rotation curves of galaxies. Potential- 
density pairs (e.g., Hernquist 1990, Dehnen 1993, Zhao 1995) have long been acknowledged to be very useful in model building 
and in checking numerical simulations. It is even more interesting to find simple spherical and non-spherical models in MOND 
and TeVeS. We explore a spherical model with a scalar field of the form 

0W-o 2 m(i + ^), (ii) 

where p is a dimensionless number, and 

,,2 



^0 2 



r = — , v- = VGM a . (12) 
a 

Then we have, according to the equation for the scalar field (Eq. |S), that 

M [r„(r + r„)-rg/4] 1 

P(r) = jj-2, r h = (p+ -)r (13) 

2nr [(r + r h ) 2 - rjj/4] 2 

The density profile is shown for p=0.5 and 2.0 in Fig.l. This MONDian density distribution is realistic: the model mimics the 
Sersic profile of an elliptical galaxy. It also mimics a Hernquist (1990) model with scalelength r^: 

Ph{r) = - — r — rg- (14) 

Airr(r + rh) 

We can then also easily derive g s = Vo/(r + ro +pro), and from Eq.© calculate the Newtonian and total potential of the 
model: 

$jv(r) = $-0(r), $(r)=«gln( 1 + — ) . (15) 

\ pro J 

Thus, for the gravity we have 

9(r) = ^ = (16) 
v ' r r + pro 

The corresponding matter density (baryons + DM) in the Newtonian framework is 
M„(l + 3Ep) 

P (0 + ^M(r)= 47rro(pro + r)2 . (17) 

v 2 

The asymptotic circular velocity is vo, and the maximum gravity go occurs at small r (near the origin), where g(r) — = ^ 1 
(note that in the flattened models of §7 go will be related to p by go = —). Given this, it is interesting to show the spherical 
model with p=0.5 and 2.0, hence the intermediate gravity range go = 2ao and ao/2 (Figs. 1-4). Most bulges and ellipticals, 
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Figure 1. Shows the density of Eg. 1131 for p=0.5 and 2.0. tq and are taken to be unity, and G=0. 00442 in all following plots unless 
stated otherwise. Overplotted is a Hernquist density profile (cf. Ea. 1141 for p=0.5 and 2.0. It shadows the curve for p=2 from our density 
profile extremely well, but there is a slight deviation between the p=0.5 curves. 



gravitational lenses and galaxy clusters should go through this intermediate regime, while dwarf spheroidal galaxies have very 
low go/aa = ^. Although a detailed fit to observed light profile is beyond the scope of the present paper, we will show an 
application in galaxy clusters in §5. 



4.1 The baryonic distribution function 



It is interesting to ask if the MONDian potential-density pairs presented in the previous section can be realised by some 
equilibrium configurations described by certain baryonic distribution functions self-consistently. This can be done in exactly 
the same way as in Newtonian gravity. This is particularly simple for a system with constant radial anisotropy such that the 



radial velocity dispersion oy is related to the tangential dispersions by oy 
distribution F for the distribution function must take the form: 



2ol 



F(E,L) = 



ME) 



A{E) 



d(rp) 



4ttL 

Knowing from Eg. 1151 that 
r($) = pro 
we obtain 

Mo 



(/<!> 



2<t j,. In this case the baryonic phase-space 



(18) 



* 1 
exp — — 1 



F(E,L) = 



A d 2A 2 + i 



1)A + 3p(2p + 1)) 



16ir 2 r%v%L 



p(p + A) 3 



- E 

A = exp — 7t- . 



(19) 



(20) 



For this radially anisotropic model we can calculate the radial velocity dispersion by solving the Jeans equation as done 
for Newtonian system (e.g., Angus & Zhao 2006) 



2a 2 e 



2ai 



rp 



Vcpdr. 



Substitute in the expression for the circular velocity V c we have 



al{r) 



i 



4(2p+ l)s-4 



(f + !)(.»: 



+ 



S + 1 S + 1 



+ ■ 



+ 



is- 



(21) 



(22) 



where s = l + 2p+ A likewise expression can be found by calculating the moments from the distribution function F(E, L). 

An interesting feature of this model is that it has nearly isothermal velocity dispersions everywhere for any value of the 
parameter p. Fig.(|5J shows oy(r) for the two limiting cases of p = and p = oo. To understand this note that at large 
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p = 



p = inf iriit y 



Figure 2. Shows the radial velocity dispersion of anistropic model given by Eo. 1221 for the limiting cases p=0 and p=oo. Here no and 
7*0 are unity. 




Figure 3. Shows the anisotropic distribution function from Ea, l20l at L = 1 for p=0.5 and 2.0 

radii p ~ r~ 4 , V c — > Wo and of — > At small radii, rp — > est, and (7^(0) = (p + 1)vq [3/4 + p/2 + p(p/2 + 1) In = 
0.75^, 0.470?Jo,0.421wg,0.383«o,0.3337jg for p = 0,0.5, 1,2, oo. So for all plausible values of p, of is very comparable at the 
center and at large radii. The anisotropic distribution function as a function of E is presented in Fig.Q for p = 0.5 and p= 
2.0. We can also numerically integrate the isotropic distribution function via Eddington's equation (see appendix and Fig0J. 



4.2 Gravitational lensing 

Light bending in TeVeS works very much like in General Relativity. For a ray of impact parameter R from the spherical lens, 
the bending angle 9 is given (cf. Zhao, Bacon, Taylor, Home 2006) by the following integration along the line of sight 



e(R) = 



2g± dz 



R 

g±(R,z) = g(r) — , 



(23) 



where g±(R,z) is the gravity perpendicular to the line of sight, and g(r) is the centripetal gravity at radius r = \/R 2 + z 2 . 
Using the expression for g(r) from our model (Eq. llOL we find that the deflection angle is given by 



: arctan . 



2^- when y = > 1 
y+l " pro 



e(R) 



when y 



-PI tanh-\ 

V y +1 



R 

pro 

— when y = -S. 

t> pro 



= 1 
< 1 



(24) 
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Figure 4. Shows the isotropic distribution function derived in the appendix for p=0.5 and 2.0. 
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Figure 5. Shows the bending angle in arcseconds as a function of impact parameter R and the circular speed V c in 10 2 kms 1 for 
pro=0.5Mpc and vq=949 km s^ 1 . 



FigQ>] shows the predicted bending angle as function of the impact parameter R; cases are shown for p — 0.5 and 2.0. The 
bending angle increases with the impact parameter and starts to level off beyond R = pro- This is consistent with the 
expectation of a flat rotation curve at large radii (Fig. Given the distance to the lens Di, the distance to the source D„ 
and the lens-source distance Di s we can define an effective distance as D e g — DiDi s /D„. Using this we can compute the 
convergence k and the tangential shear 7 from 

k(R) = ^- 7 (i*) 

R 2 - 2(pr ) 2 (fl - e{pr ))D cS 



2R 2 - 2(pr ) 2 R 



and they are plotted for a toy cluster in Fig.JSJ. 
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Figure 6. Shows the convergence re and tangential shear 7 for a toy cluster at an effective distance D c g=400Mpc with i>o=949kms~ 
and pro=0.5Mpc. 
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Figure 7. Shows contours of a double-centred (panel a) and a triple-centred (panel b) analytical model for the potential <E> (dashed red 
lower) and projected lensing convergence re (thick blue shaded lower) of a system resembling the bullet cluster (Clowe et al. 2004). Also 
shown are contours (spaced in 0.5 dex, in internals of factor of two) the matter volume density in Einstein-Newton gravity (upper thin 
black) and baryonic matter volume density in the classical MOND approximation (upper shaded zones). Note a density minimum for 
'MONDian' baryons at (0,0) in panel (a), and a density maximum in panel (b) where the baryons are primarily in a disky component in 
the middle. 



5 MULTI-CENTRED POTENTIALS AND THE MERGING BULLET CLUSTER 

An obvious way of obtaining flattened models is to consider a double-centred potential, whose axis of symmetry is perpendicular 
to the line of sight, i.e., a merging system viewed edge-on. In a cartesian coordinate system (x,y,z), if the two centres are 
located at (—xq, 0,0) and (a;o,0, 0), then the double-centred potential can be written as 

= ¥ + (26) 
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where n = yjx + xq) 2 + y 2 + z 2 , r 2 = y/ (x — xq) 2 + y' 2 + z 2 , and "!> is the spherical potential of Eq.(16). Here we have simply 
added two spherical analytical potentials. For example, to mimic a merging galaxy cluster we set p = 2, xo = pro = 0.5Mpc, 
wo = 949 km/s, and M = 5.21 x 1O 13 M . 

The gravitational field of such a model is easily computed by the simple superposition of two spherical fields, i.e., 

g = !!°_£l ^L_£l. ( 27) 

ri + pr 2v\ r 2 + pr 2r 2 

Likewise the light bending angle vector is the superposition of the bending angle vectors of two spherical models, so its x and 
y direction components are respectively 



?dcx\->-,y)— 2R t 2Rr> 
g, ( x y \- + 



where Ri = ijjx + xo) 2 + y 2 , R 2 = \J {x — Xo) 2 + y 2 , and 8 is the spherical bending angle of Eg. 12411 . From there we can 
estimate the projected lensing convergence 

, > 1 _. f dOdcx , d9 dcy \ k(Ri) k{R 2 ) , ofV , 
K(x,y) = -D eS -5 — + =-i— i + ( 29 ) 



2 I 9i <9j/ 

where k(_R) is given by Eq, 125H for spherical systems. Assuming D c a = 400Mpc, the resulting convergence contours are plotted 
on Fig. |7[a). This map resembles the convergence map derived from the weak lensing shear field around the merging bullet 
cluster IE 0657+56 (Clowe et al. 2004). 

So far the result does not differ from Einsteinian gravity. The difference is in the underlying matter distribution. In 
Einsteinian gravity, the convergence map is simply proportional to the surface density map, and the volume density of matter 
is the simple addition of two spherical models each with density given in Eg. 1171 . As we can see from Fig|7]there is a one-to-one 
relation between features in the convergence map and the features in the matter (dark plus baryons) distribution if the gravity 
is Einsteinian, i.e., what we see in lensing is what we have. 

The situation in MOND/TeVeS is different. In order to properly derive the corresponding density in TeVeS, one would 
need to know what part of the double-centred potential is due to the scalar field. This will be done in the limiting case of 
scale-free flattened models in §7. Here we use the classical MOND approximation. From the gravity or the potential we can 
then directly calculate the corresponding baryonic isodensity contours using Eq.(2) and Eq.(9) with a = 1. The MONDian 
baryonic matter has a rich structure. Fig[7{a) shows that at the centre (x = 0, y = 0) the MONDian volume density reaches a 
local minimum while the convergence map shows a saddle point. This cautions against a naive deprojection of the convergence 
map in MOND/TeVeS. 



5.1 Triple-centred baryonic system and the bullet cluster IE 0657+56 

In the case of the bullet cluster, there are three baryonic mass concentrations. Clowe et al. (2004) argue that the bulk of 
baryonic mass is in the form of X-ray gas, which is shocked and displaced from the two optical centres of the colliding binary 
cluster. It was argued that lensing in any MONDian theory should produce shear maps centred on the dominating X-ray gas 
instead of the lesser baryonic mass responsible for the optical light. The fact that the convergence map coincides with the 
two optical centres is presented as direct evidence for the presence of collisionless dark matter, unaffected by the shock, and 
respecting the optical centres. 

Of course, this is not so surprising since it is well known that MOND still needs an unseen matter component in galaxy 
clusters (Sanders 2003). But, in the case of the bullet, a key element of the line of reasoning is that the geometry of the lensing 
map in TeVeS reflects the underlying baryons even in highly non-spherical geometries. To illustrate the range of possibilities 
in triple-centred systems, let us consider the following potential 

$tc(x,y, z) = [fea + (1 - fei - k 2 )H{x)} $(n) + [fa + (1 - fa - k 2 )H{-x)} #(ra). (30) 

The terms involving the Heaviside-function create the effect of a razor thin disk at x — (reminiscent of the well-known 
Kuzmin disk), with a sudden change of gravity in the x-direction at the midplane x = 0. The potential is continuous across 
the x — plane. The deflection angle and convergence map remain those of the superposition of two spherical deflectors, 
although there is now a sudden change in the weighting of the two deflectors as one crosses the midplane x — 0. Finally one 
can apply the MONDian Poisson equation (Eq. 2) to derive the baryonic density. 

Here we choose fa = k 2 = 0.2 so that we have a prominent baryonic component in between the two cluster centres. 
Fig|7Jb) shows a rather regular looking convergence map, but a complex MONDian baryonic density distribution, somewhat 
resembling that of the bullet cluster. In particular, there is now a ridge of matter centred oni = and two irregular baryonic 
components more or less centred on the x = ±xo- This example provides evidence that a regular-looking convergence map is 
not incompatible with a MONDian multi-centred baryonic mass distribution. 

As argued in Zhao, Bacon, Taylor & Home (2006) and Zhao & Qin (2006) the convergence k can be non-zero where there 
is no projected matter in MOND/TeVeS, something that is not possible in Einsteinian gravity. This implies that a lensing 
convergence map does not simply translate a baryonic surface density map in MOND. Our models here show the range of 
non-trivial baryonic geometries for multi-centred potentials in MOND. Although the models here are of only indicative value 
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Figure 8. Shows a flattened non-scale-free potential with go = arj/2 (p=2), e=0.25, upper panel shading showing density contours (for 
a = 1), bottom dashed showing Hernquist density profile (cf. Eg, 1141 with axis ratio 2/3. 



(as the approximation of classical MOND has been used, see §7 for typical differences with TeVeS) it does caution drawing 
conclusions about TeVeS before performing a careful lensing analysis of the bullet cluster (Clowe et al. 2004) in the framework 
of TeVeS. 



6 FLATTENED POTENTIALS AND ORBITS 



The spherical potential of §4 can also be generalised into a flattened or triaxial potential. For example, we can consider a 
model with an axisymmetric potential 



<£>M) = |ln 



r \ To 

H + 2e cos 

pro J pro 



(31) 



where 9 is the angle with the z-axis: e > corresponds to an oblate potential and e < to a prolate potential. 

The corresponding density in the classical MOND of Bekenstein & Milgrom (1984) can be calculated by feeding this 
potential into the modified Poisson equation (Eq.|2J in spherical coordinates: 



47rGp(r, 9) = 



r 2 dr dr 



+ 



d fi sin #<9$ 
lindde d6 ■ 



(32) 



The expression for p can be obtained analytically, but the general expression is too tedious to be given here. Nevertheless we 
can calculate orbits in this potential numerically. The stars in this flattened MOND potential are typically on loop orbits as 
in a flattened Newtonian potential (Fig. 0. We shall concentrate hereafter on a limiting case of this flattened model, in order 
to compare the predictions of the classical MOND gravity and of TeVeS. 



7 SCALE-FREE FLATTENED MODELS: MULTI-FIELD THEORY VS. ONE-FIELD THEORY 

In the limit r/pro — > 0, the previous potential reduces to a scale-free form when expanding Eq. l3H . We shall now consider 
such models with a total potential 

3>(r, ff) = rg r (cos 9), where g r (c) = (1 + ec 2 )go, (33) 
and for convenience we define 

c = cos9. (34) 

Clearly, go = vg/pro is thus the gravity along the major axis in the c = (or 9 — 7r/2) plane, and the parameter e leads to 
flattening (e = reduces to the spherical case). The parameter go is NOT the acceleration constant ao of MOND, but is a 
parameter of the model linked to it through the definitions (see Eq. I12H of vo and ro (go = ao /p) ■ 
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This class of scale-free models corresponds to flattened galactic bulges, ellipticals, dwarf spheroidal galaxies, or centre 
of galaxy clusters depending on the value of the equatorial gravity go- The radius-independent gravitational force in these 
models is: 

S (c)=V<S>(r,9) = g r (c)r + g g (c)§, (35) 

2,1/2 ( dgr\ _ , 2x1/2 



where g g {c) = (1 - c 2 ) 1/2 [^) = 2ce 5o (l - c ) I/2 . (36) 
The amplitude of gravity g is defined as 

g{c) = [g r {cf +ge{c) 2 ] 112 . (37) 

We are now going to explore the density corresponding to this potential in four different frameworks: (i) Newtonian 
gravity where the computed density corresponds to the baryonic and DM distribution, (ii) classical one-field MOND gravity 
(Bekenstein & Milgrom 1984, approximation used in §5), (iii) TeVeS multi-field theory, (iv) original MOND formula (Eq. 
approximation. 1 

For the first time in the literature, we provide here a quantitative comparison of the three ways of describing the 
MOND phenomenology without invoking DM (and of Newtonian gravity with DM), by comparing the underlying density 
corresponding to the potential of Eg. 1331 , 

7.1 Newtonian gravity with dark matter 

Assuming Newtonian gravity, the density (corresponding to the stellar and dark densities) is 

p(r,6) = (47rG) _:L V 2 $(r,0) = (AnGr)- 1 L(c), (38) 

where L(c) = 2g r (c) + -^[(1 - c 2 )g' r (c)] = C[g r (c)\ = 2(1 + e)g - 4e 5o c 2 , (39) 

We have defined the differential operator £ = 2+ ^(1 + c 2 )^. This toy model has a r _1 cusp, a rising rotation curve V oc ■s/r, 
and M(r) oc r 2 (similar to uniform disks). 

7.2 One-field MOND gravity 

Now, assuming a classical MOND gravity (as we did in §5) with an interpolating function jJ,(g), the underlying (purely 
baryonic) density can be computed from 

PmM) = (4 7 rG)~ 1 V.(MV$(r,e)) = (^Gr) -1 L M (c), (40) 
where L M (c) = p(g)L(c) + (1 - c 2 )^M^l, (41) 

where L(c) is given by Ea. l|39^ and fi(g) by Eq.0 with g as in Eq. jjffi . This baryonic density pM(r,8) is compared with 
the baryonic+DM density of a Newtonian model with the same potential in Fig|5J for different values of go and different 
interpolating functions. 

7.3 Multi-field TeVeS gravity 

To compute the baryonic density in §5 when making a toy model of the bullet cluster, we used the one-field MOND gravity. 
We are now going to show the differences that could be expected when using the multi-field TeVeS instead. The problem in 
that case is that we must find the relative contribution of the scalar field and of the Newtonian potential to the total potential 
before computing the underlying baryonic density. 

Assuming a multi- field gravity (see Eqs. 4 and 5) with a scalar field cj> and a Newtonian potential $jv = $ — <j>, we have 
that the scalar field can be written in the form 

4>{r,9) = rg 3ir (cos0), where g 3 , r (c) = g r (c) - g n ,r{c) = g (l + ec 2 ) - g n , T {c), (42) 

where g n ,r and g a , r are respectively the Newtonian gravity in the r direction and the scalar field strength in the f direction, 
both generated by the common density distribution we are looking for, p s (r,6). 
The total scalar field strength is 

g.(c) = \V<P(r,6)\ = (g 2 sAc)+gle(c)) 1/2 , where g s ,e = (1 - c 2 ) 1 ' 2 (fe) . (43) 



1 In a flattened system, Eq.(l) does not derive from a proper theory of gravity (energy is generally not conserved). The only way to 
precisely recover it is to consider purely circular orbits in axisymmetric systems for non-relativistic modified inertia toy models (see 
Milgrom 1994). In a flattened system, the corresponding formula for the classical one-field MOND gravity will differ from Eq.(l) by a 
curl-field term (Eq. while this curl-field will only affect the scalar field in TeVeS. 
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Table 1. Shows the correction to the rotation curve due to the curl field for six o's (cf. Eq. and for three representative gravity 
strengths (go = Qo/2, ao and 2ao). e = 0.8. 
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-0.274 
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0.002 


0.002 


0.000 


0.007 


-1.238 


1.244 


-0.731 


-7.074 


6.343 


1.0 


-0.114 


-1.985 


1.871 


-0.255 


-0.102 


-0.152 


-0.820 


-2.132 


1.311 



The baryon density is then given by both the scalar and Newtonian Poisson equations: 

p s {r,6) = (AnGry'MgsAc)] + (1 - c 2 ) ^l^ }, (44) 
and 

Ps (r,0) = {A-kGt)- 1 C[g r {c) - g s , r (c)] (45) 

Combining the above two equations, we use the numerical relaxation method to solve the following second-order ODE for 
fl«,r(c): 

[1 + p s ]£[g s Ac)] + (1 - c 2 ) dfls d [ ° s) ^ = He). (46) 

We set the boundary conditions such that the solution is regular at c = ±1. From this we derive the baryonic density p s (r, 6). 

We plot the equal density contours predicted from the three theories for a model with go — ao/2 and 2ao, a=0.0 and 1.0 
and e=0.8 in Fig|5| Note that some of the models have an unphysical heart-shaped density distribution. This doesn't prevent 
us comparing them in order to show the differences between the different formulations of MONDian modified gravity theories. 



7.4 MOND formula approximation 

In the case of the multi-field TeVeS gravity, if one subtracts the true gravity go in the equatorial plane from the one that we 
would derive when applying Eq.(l) to the Newton-Einstein gravity, we obtain the correction to the rotation curve due to the 
curl field: 

9s AO) (l + T^£w)-5o 

A 3 = i ——1 . (47) 

go 

When this number is negative, the true MOND force is larger than the one derived using Eq.(l), while the opposite is true 
when this number is positive. 

In the case of the one-field MOND gravity, the actual Newtonian gravitational force <7iv(c) can be calculated by solving 
the Newtonian Poisson equation for the density pM(r,9). Then, replacing <? s ,r(0) by go — <?iv(0) in Eq. l47H yields Am, the 
correction to the rotation curve due to the curl field in the non-relativistic MOND. Although the two values do not correspond 
to the same mass density, the relatively large values of the difference imply non-negligible differences. Those values are listed in 
Table 1 for three different values of go (ao/2, ao, 2ao), corresponding to the intermediate gravity regime, and for many different 
interpolating functions. We thus showed that it is necessary to do rigourous calculations using the scalar field prescription if 
one wants to address the intermediate regime quantitatively in complex geometries. 



8 CONCLUSION 

In summary, we have presented analytical models of dynamics and lensing in MOND/TeVeS. 

1. We proposed (iJSJ a useful set of interpolating functions for MOND, with a physical counterpart in TeVeS, con- 
trary to the standard interpolating function commonly used to fit galactic rotation curves. 

2. Using our interpolating functions, we found a useful family of spherical (i0 models in MOND/TeVeS, with mod- 
erate 1/r cusps. Those models were then extended to non-spherical ones, by considering multi-centred models flattened 
models (0 and scale-free flattened models (50. 
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Figure 9. Shows isopotcntials (dotted line) for the potential of Ea . 1331 . for 2 values of the equatorial gravity go (go = &0 /2 and go = 2ao), 
and 2 values of a (see Eq.(9)). The full black line corresponds to isodensity contours in Newtonian gravity (baryons+DM). The dashed 
line corresponds to baryonic isodensities for the one-field MOND and the dot-dashed line to the baryonic isodensities for the multi-field 
TeVeS. For the same potential, the typical density in TeVeS is slightly higher than in MOND near the centre, and slightly lower in the 
outskirts. 



3. We showed that the lensing and the orbits in these spherical or flattened models are rather similar to the expec- 
tation in Einstein-Newton gravity, but still trigger a few surprises in extreme geometries. In multi-centred models, the 
convergence map does not always reflect the projected matter in the lens plane in MOND. This cautions simple interpretations 
of the analysis of weak lensing in the bullet cluster IE 0657-56 (Clowe et al. 2004, see Fig. 0. 

4. In flattened scale-free models we also found that there are differences in the potential-density pairs between 
the single-field MOND formulation of Bekenstein & Milgrom (1984) and the multi-field TeVeS formulation (Bekenstein 2004). 
However, the differences are probably not sufficient to solve the MONDian mass discrepancy in galaxy clusters (Sanders 
2003) as suggested by Bekenstein (2005), unless other parameters of the TeVeS theory, e.g., S (see Eq.|lJ, are important. 
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Note that our results are not inconsistent with the well-known result that rotation curves of astronomical disks are 
insensitive to the detailed formulation of MOND. While the very squashed axisymmetric systems of J3are not very realistic 
by themselves, they seem to suggest that the three versions of MOND are likely to show greater differences in systems 
of complex multi-centred geometry, which is realistic for systems undergoing mergers. 2 Altogether we would argue that 
the most promising systems to test different versions of MOND are systems of lower symmetry for which MOND was not 
designed and the least is known. We have shown here that application of this test requires highly non-trivial computation to 
be done properly. 
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APPENDIX A: ISOTROPIC DISTRIBUTION FUNCTION 



Here we present the full isotropic distribution function corresponding to the models of ^ We start from the Eddington's 
formula (see Binney & Tremaine 1987) 



F(E) 



V8tt 2 



d 2 p d$ 



+ 



, , d$ 2 v$ - E v$~e Vd$;*= ( 

Initially we define p as a function of $ 

_ M A 2 ((2p+l)A- 1 + l) _ -£ 

P( ) 47rp 2 r3(A- 1 -l)(pA- 1 + l) 2 ' 6XP v 2 Q ' 

The first derivative of p w.r.t. $ is 

2p(2p + 1)A" 1 + (3p 2 - 3p - 1) + (3p - 1)A + A 2 



27rro«o dp 
M d$ 



(A- 1 - l) 2 (pA"! + l) 3 



(Al) 



(A2) 



(A3) 



Evaluated at $ = oo we get ^|$ =00 = 0. The second derivative w.r.t. <!> is 



2 (2p + 1)A -3 + p(-23p 2 + 15p + 7)A~ 2 + (9p 3 - 29p 2 + lip + 2^^ + (12p 2 - 20p + 3) + (8p - 5)A + 2A 2 



M cM? 2 
The reduced Eddington's formula 



P 2 (A- 



l)3(pA-i + l) 4 



(A4) 



2 For example, satellites in the tidal field of a galaxy have interesting Roche lobe shapes, which contain detailed information on the 
different laws of gravity (Zhao & Tian, 2006). 
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F{E) 



d 2 p d$ 



d$ 2 V* - E 



(A5) 



can now easily be numerically integrated as shown in Fig.lQJ. This paper has been typeset from a TgX/ IMpjX file prepared 



by the author. 



